output_version <- '2024-07-22'

## fishnet
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"

#For each cell (same cell, 50km, 100km), we need the area shares of the various categories below:
require(exactextractr)

sum_cover <- function(x){
  list(x %>%
         group_by(value) %>%
         summarize(total_area = sum(coverage_area)) %>%
         mutate(sh_landsuit = total_area/sum(total_area)))
  
}

## woody biomass https://daac.ornl.gov/daacdata/global_vegetation/Gridded_Biomass_Africa/
## Hanan, N.P., L. Prihodko, C.W. Ross, G. Bucini, and A.T. Tredennick. 2020
woody_biomass_local_fn <- file.path(wkdir,'local/gis_data/007_environment/woody_biomass/woody_biomass.tif')
xfun::dir_create(dirname(woody_biomass_local_fn))

if(!file.exists(woody_biomass_local_fn)){
  woody_biomass_fn <- file.path(wkdir,'woody_biomass.tif')
  woodbiomass_prj <- terra::rast(woody_biomass_fn)
  
  # project geographic 1km
  extent_geo <- ext(c(-180, 180, -90, 90))
  r1km <- rast(extent_geo, res=0.0083333333)
  woodbiomass_global <- terra::project(woodbiomass_prj, r1km, method="near", threads = 16)
  crs(woodbiomass_global) <- "+proj=longlat +datum=WGS84"
  fishnet_ext <- c(0,25,0,24)
  woodbiomass <- crop(woodbiomass_global,fishnet_ext)  
  writeRaster(woodbiomass,woody_biomass_local_fn)
} else {
  woodbiomass <- rast(woody_biomass_local_fn)
}

##
## RADIUS BUFFER 0 (grid cell only) 50KM (FROM POINT), 100KM (FROM POINT)
##

## memory issues
n_objectid_list=unique(fishnet_adm0_only$OBJECTID)
total_chunks <- split(n_objectid_list, ceiling(seq_along(n_objectid_list)/1000))
length(total_chunks)

  ## memory issues
  require(doParallel)
  
  for(chunk_index in 1:length(total_chunks)) {
    print(chunk_index)
    woodybiomass_chunk_out_fn <- file.path(wkdir,'proc_data',
                                       'woody_biomass_chunks',sprintf('WOODYBIOMASS_%s_of_%s_%s.csv',chunk_index,length(total_chunks),output_version))
    xfun::dir_create(dirname(woodybiomass_chunk_out_fn))

    if(!file.exists(woodybiomass_chunk_out_fn)){
      n_objectid_list_chunk <- as.numeric(unlist(total_chunks[chunk_index]))
      
      n_cores <- min(10,length(n_objectid_list_chunk))
  
      cl <- makeCluster(n_cores)
      registerDoParallel(cl)
      
      results <- foreach(objectid_j=n_objectid_list_chunk,
                         .export=c(ls()),
                         .packages = c("dplyr", "sf","exactextractr","spatialEco","terra")) %dopar% {
                           # Modulus operation
                           if(objectid_j %% 100==0) {
                             # Print on the screen some message
                             cat(paste0("iteration: ", objectid_j, "\n"))
                           }
                           
                           fishnet_lake_chad_ntl_adm0_x_d01_shp_fn <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01.shp')
                           fishnet_adm0_only_j = read_sf(fishnet_lake_chad_ntl_adm0_x_d01_shp_fn, query = sprintf('SELECT * from fishnet_lake_chad_ntl_adm0_x_d01 WHERE OBJECTID = %s',objectid_j))
                           st_crs(fishnet_adm0_only_j) <- "EPSG:4326"
                           
                           if(dim(fishnet_adm0_only_j)[1]>0){
                             
                             fishnet_adm0_only_buf100km_bbox <- round(st_bbox(spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=150000)),0)
                             ## crop raster to 100km buffer
                             woody_raw_i <- terra::rast(woody_biomass_local_fn)
                             
                             woody_crop_j <- terra::crop(woody_raw_i,fishnet_adm0_only_buf100km_bbox,snap="out")
                             all_na <- all(is.na(values(woody_crop_j)))
                             if(all_na==FALSE){
                               ## 0km
                               x <- exact_extract(woody_crop_j, fishnet_adm0_only_j, 'mean')
                               
                               #merge the list elements into a df
                               woodybiomass_0 <- as.data.frame(cbind(x, objectid = objectid_j)) %>%
                                 mutate(buffer_km = 0)
                               
                               rm(x)
                               # projected 50km
                               fishnet_adm0_only_buf50km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=50000)
                               x <- exact_extract(woody_crop_j, fishnet_adm0_only_buf50km, fun = 'mean')
                               
                               #merge the list elements into a df
                               woodybiomass_50 <- as.data.frame(cbind(x, objectid = objectid_j)) %>%
                                 mutate(buffer_km = 50)
                               
                               # projected 50km
                               fishnet_adm0_only_buf100km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=100000)
                               x <- exact_extract(woody_crop_j, fishnet_adm0_only_buf100km, fun = "mean")
                               #merge the list elements into a df
                               woodybiomass_100 <- as.data.frame(cbind(x, objectid = objectid_j)) %>%
                                 mutate(buffer_km = 100)
                               
                               out <- rbind(woodybiomass_0,woodybiomass_50,woodybiomass_100)
                               names(out)[1] <- "woody_biomass"
                               
                               # add objectid
                               out$objectid <- objectid_j
                               
                                 
                               # clean-up
                               rm(x,fishnet_adm0_only_buf100km,fishnet_adm0_only_buf50km)
                               
                               return(out)  
                             }
                           }
                           
                         }
      
      parallel::stopCluster(cl)
      
      system.time(woodybiomas_dt <- data.table::rbindlist(results))
      woodybiomas_df <- woodybiomas_dt %>% as.data.frame() %>%
          dplyr::select(objectid,woody_biomass,buffer_km) %>%
          arrange(objectid,buffer_km)
      
      data.table::fwrite(woodybiomas_df,woodybiomass_chunk_out_fn)
      
      # clean-up
      rm(woodybiomas_df)
      rm(woodybiomas_dt)
    }
      
    
  }
  
  # merge all chuncks
  woodybiomass_df_files <- sort(list.files(dirname(woodybiomass_chunk_out_fn),
                                        pattern=sprintf("WOODYBIOMASS"),
                                        full.names=T))
  tables <- lapply(woodybiomass_df_files, function(z) read.csv(z))
  woodybiomass_df_all <- data.table::rbindlist(tables)
  
  #
  woodybiomass_out <- woodybiomass_df_all %>%
    dplyr::select(objectid,woody_biomass,buffer_km) %>%
    arrange(objectid,buffer_km)
  
  # out file  
  woodybiomass_out_fn <- file.path(wkdir,'proc_data',sprintf('WOODYBIOMASS_%s.dta',output_version))
  
  # Adding the label value
  attr(woodybiomass_out$objectid, "label") <- "objectid of grid"
  attr(woodybiomass_out$woody_biomass, "label") <- 'Mean Woody Biomass 1-km resolution in megagrams per hectare (Mg ha-1) : Hanan et al.'
  attr(woodybiomass_out$buffer_km, "label") <- 'Buffer'
  attr(woodybiomass_out, "label") <- sprintf("woody_biomass.R last run on %s",Sys.Date())
  head(woodybiomass_out)
  summary(woodybiomass_out)
  
  ## Save the data
  haven::write_dta(woodybiomass_out, woodybiomass_out_fn)
  